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Abstract — We consider the problem of exploiting the mi- 
crogenerators dispersed in the power distribution network in 
order to provide distributed reactive power compensation for 
power losses minimization and voltage support. The proposed 
strategy requires that all the intelligent agents, located at the 
microgenerator buses, measure their voltage and share these data 
with the other agents on a cyber layer. The agents then actuate the 
physical layer by adjusting the amount of reactive power injected 
into the grid, according to a feedback control law that descends 
from duality-based methods applied to the optimal reactive power 
flow problem subject to voltage constraints. Convergence of the 
algorithm to the configuration of minimum losses and feasible 
voltages is proved analytically for both a synchronous and an 
asynchronous version of the algorithm, where agents update their 
state independently one from the other. Simulations are provided 
in order to illustrate the algorithm behavior, and the innovative 
feedback nature of such strategy is discussed. 

I. Introduction 

Recent technological advances, together with environmental 
and economic challenges, have been motivating the deploy- 
ment of small power generators in the low voltage and medium 
voltage power distribution grid. The availability of a large 
number of these generators in the distribution grid can yield 
relevant benefits to the network operation, which go beyond 
the availability of clean, inexpensive electrical power. They 
can be used to provide a number of ancillary services that are 
of great interest for the management of the grid [1], [2]. 

We focus in particular on the problem of optimal reactive 
power compensation for power losses minimization and volt- 
age support. In order to properly command the operation of 
these devices, the distribution network operator is required 
to solve an optimal reactive power flow (ORPF) problem. 
Powerful solvers have been designed for the ORPF problem, 
and advanced optimization techniques have been recently 
specialized for this task [3], [4]. However, these solvers assume 
that an accurate model of the grid is available, that all the grid 
buses are monitored, that loads announce their demand profiles 
in advance, and that generators and actuators can be dispatched 
on a day-ahead, hour-ahead, and real-time basis. For this 
reason, these solvers are in general offline and centralized, and 
they collect all the necessary field data, compute the optimal 
configuration, and dispatch the reactive power production at 
the generators. 

These tools cannot be applied directly to the ORPF problem 
faced in microgrids and, more in general, in low/medium 
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voltage power distribution networks. The main reasons are 
that not all the buses of the grid are monitored, individual 
loads are unUkely to announce they demand profile in advance, 
the availability of small size generators is hard to predict 
(being often correlated with the availability of renewable en- 
ergy sources). Moreover, the grid parameters, and sometimes 
even the topology of the grid, are only partially known, and 
generators are expected to connect and disconnect, requiring 
an automatic reconfiguration of the grid control infrastructure 
(the so called plug and play approach). 

Different strategies have been recently proposed in order to 
address these issues, ranging from purely local algorithms, 
in which each generator is operated according to its own 
measurements [2], to distributed approaches that do not require 
any central controller, but still require measurements at all the 
buses of the distribution grid [5]. Only recently, algorithms 
that are truly scalable in the number of generators and do not 
require the monitoring of all the buses of the grid, have been 
proposed for the problem of power loss minimization (with 
no voltage constraints) [6], [7]. While these algorithms have 
been designed by specializing classical nonlinear optimization 
algorithms to the ORPF problem, they can also be considered 
as feedback control strategies. Indeed, the key feature of these 
algorithms is that they require the alternation of measurement 
and actuation based on the measured data, and therefore they 
are inherently online algorithms. In particular, the reactive 
power injection of the generators is adjusted by these algo- 
rithms based on the phasorial voltage measurements that are 
performed at the buses where the generators are connected. 
The resulting closed loop system features a tight dynamic 
interconnection of the physical layer (the grid, the genera- 
tors, the loads) with the cyber layer (where communication, 
computation, and decision happen). In this paper, we design 
a distributed feedback algorithm for the ORPF problem with 
voltage constraints, in which the goal is to minimize reactive 
power flows while ensuring that the voltage magnitude across 
the network is larger than a given threshold. 

In Section III, a model for the cyber-physical system of a 
smart power distribution grid is provided. In Section IV, the 
optimal reactive power flow problem with voltage constraints 
is stated. An algorithm for its solution is derived in Section V, 
by using the tools of dual decomposition. A synchronous and 
an asynchronous version of the algorithm are presented in 
Section VI and Section VII, respectively. The convergence 
of both the proposed algorithms is studied in Section VIII. 
Some simulations are provided in Section IX, while Section X 
concludes the paper discussing some relevant features of the 
feedback nature of the proposed strategy. 
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II. Mathematical preliminaries and notation 

Let Q = (V,f,(T, r) be a directed graph, where V is the 
set of nodes, £ is the set of edges, and ct, t : £ — ^ V are two 
functions such that edge e E £ goes from the source node 
a{e) to the terminal node r(e). 

Given two nodes of the graph h,k e V, we define the 
path Vfik — {vi, ■ ■ ■ ,ve) as the sequence of nodes, without 
repetitions, such that 

• vi = h 

• Vi = k 

• for each i = 1, ...,£— 1, the nodes Vi and u^+i are 
connected by an edge. 

In the rest of the paper we will often introduce complex- 
valued functions defined on the nodes and on the edges. These 
functions will also be intended as vectors in C" (where n ~ 
|V|) and C'^l. Given a vector u, we denote by u its (element- 
wise) complex conjugate, and by its transpose. We denote 
by 9fi(w) and by the real and the imaginary part of u, 

respectively. 

Let moreover A e {0,±1}'^'^" be the incidence matrix of 
the graph Q, defined via its elements 

{—1 if V — a{e) 
1 ifv = T{e) 
otherwise. 

If the graph Q is connected (i.e. for every pair of nodes 
there is a path connecting them), then 1 is the only vector in 
the null space ker A, 1 being the column vector of all ones. 
We define by l^, the vector whose value is 1 in position v, 
and everywhere else. 

III. Cyber-physical model of a smart power 

DISTRIBUTION GRID 

In this work, we envision a smart power distribution net- 
work as a cyber-physical system, in which 

• the physical layer consists of the power distribution 
infrastructure, including power lines, loads, microgenera- 
tors, and the point of connection to the transmission grid, 
while 

• the cyber layer consists of intelligent agents, dispersed 
in the grid, and provided with actuation, sensing, com- 
munication, and computational capabilities. 

A. Physical layer 

For the purpose of this paper, we model the physical layer 
of a smart power distribution network as a directed graph Q, 
in which edges represent the power lines, and nodes represent 
the devices that are connected to the grid (see Figure 1, middle 
panel). These include loads, microgenerators, and also the 
point of connection to the transmission grid (called point of 
common coupling, or PCC, and indexed as node 0). 

We limit our study to the steady state behavior of the system, 
when all voltages and currents are sinusoidal signals at the 
same frequency. Each signal can therefore be represented via 
a complex number y ~ \y\eJ^y whose absolute value \y\ 
corresponds to the signal root-mean-square value, and whose 




Figure I. Schematic representation of the microgrid model. In the lower 
panel the physical layer is represented via a circuit representation, where black 
diamonds are microgenerators, white diamonds are loads, and the left-most 
element of the circuit represents the PCC. The middle panel illustrates the 
adopted graph representation for the same grid. Circled nodes represent both 
microgenerators and the PCC. The upper panel represents the cyber layer, 
where agents (i.e. microgenerator nodes and the PCC) are also connected via 
some communication infrastructure. 



phase Zy corresponds to the phase of the signal with respect 
to an arbitrary global reference. 

In this notation, the steady state of the grid is described by 
the following system variables (see Figure 1, lower panel): 

• li e C", where Uy is the grid voltage at node v; 

» i £ C", where iy is the current injected by node v; 

• ^ G Cl'^l, where is the current flowing on the edge e. 

For every edge e of the graph, we define by Zg the 
impedance of the corresponding power line. We assume the 
following. 

Assumption 1. All power lines in the grid have the same 
inductance/resistance ratio, i.e. 

Zg — 6^ I Zg I 

for any e in £ and for a fixed 9. 

This assumption is satisfied when the grid is relatively 
homogeneous, and is reasonable in most practical cases (see 
for example the IEEE standard testbeds [8]). 

The following physical constraints are satisfied by m, i and 
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communication 



A^i + ?: = 0, 

Au + e'^Z£, = 0, 



(1) 
(2) 



where A is the incidence matrix of Q, and Z = diagdzel, e S 
8) is the diagonal matrix of line impedances. Equation (1) 
corresponds to Kirchhoff's current law (KCL) at the nodes, 
while (2) describes the voltage drop on the edges of the graph. 
From (1) and (2) we can also obtain 



(3) 



where L is the weighted Laplacian of the graph L := 

A^Z-^A. 

Each node v of the grid is also characterized by a law 
relating its injected current iy with its voltage m„. We model 
the PCC as an ideal sinusoidal voltage generator at the 
microgrid nominal voltage Un with arbitrary, but fixed, angle 



J4> 



(4) 



In the power system analysis terminology, node is then a 
slack bus with fixed voltage magnitude and angle. 

We model loads and microgenerators (that is, every node v 
of the microgrid except the PCC) via the following law relating 
the voltage Uy and the current 

Uyiy = s^,, Vu e V\{0}, (5) 

where Sy is the injected complex power. The quantities 

Py := 3?(s„) and 3(s„) 

are denoted as active and reactive power, respectively. The 
complex powers Sy corresponding to grid loads are such that 
{py < 0}, meaning that positive active power is supplied to 
the devices. The complex powers corresponding to microgen- 
erators, on the other hand, are such that {py > 0}, as positive 
active power is injected into the grid. In the power system 
analysis terminology, all nodes but the PCC are being modeled 
as constant power or P-Q buses. Microgenerators fit in this 
model, as they generally are commanded via a complex power 
reference and they can inject it independently from the voltage 
at their point of connection [9], [10]. 



B. Cyber layer 

We assume that every microgenerator, and also the PCC, 
correspond to an agent in the cyber layer (see the upper panel 
of Figure 1). We denote by C (with \C\ — m) this subset of 
the nodes of Q. 

Each agent is provided with some computational capability, 
and with some sensing capability, in the form of a phasor 
measurement unit (i.e. a sensor that can measure voltage 
amplitude and angle [11]). Agents that corresponds to a 
microgenerator can also actuate the system, by commanding 
the amount of reactive power injected by that microgenerator 
(see Figure 2). 

Finally, agents can communicate, via some communication 
channel that could possibly be the same power lines (via power 
line communication - PLC - technology). Motivated by this 
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Figure 2. A schematic representation of the agents' capabilities and of the 
way in which agents of the cyber layer interface with the physical layers. The 
first panel represent the agent con'esponding to the PCC (node 0), which is 
provided with voltage measurement capabilities. The second panel represents 
all the other agents, corresponding to microgenerator nodes, which are instead 
provided with both measurement and actuation capabilities. All the agent can 
also communicate via some communication channel, and have some on-board 
computational power. 



possibility, we define the neighbors in the cyber layer in the 
following way. 

Definition 1 (Neighbors in the cyber layer). Let h ^ C be an 

agent of the cyber layer. The set of agents that are neighbors 
of h in the cyber layer, denoted as Af{h), is the subset of C 
defined as 

M{h) = {fc e C I "iVhk, VhkDC^ {h, k}} . 



Figure 3 gives an example of such set. 

We assume that every agent h E C knows its set of 
neighbors Af{h), and can communicate with them. Notice 
that this architecture can be constructed by each agent in a 
distributed way, for example by exploiting the PLC channel 
(as suggested for example in [12]). This allows also a plug- 
and-play reconfiguration of such architecture when new agents 
are connected to the grid. 

IV. Optimal reactive power flow problem 

We consider the problem of commanding the reactive power 
injection of the microgenerators in order to minimize power 
distribution losses on the power lines and to guarantee that 
the voltage magnitude stays above a given threshold. The 
decision variables (or, equivalently, the inputs of the system) 
are therefore the reactive power commands qh, h £ C\{0}. 
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Figure 3. An example of neighbor agents in the cyber layer. Circled nodes 
(both gray and black) are agents (nodes in C). Nodes circled in black belong 
to the set M{h) C C. Node circled in gray are agents which do not belong to 
the set of neighbors of h. For each agent k e JV(h), the path that connects 
htok does not include any other agent besides h and k themselves. 



Power distribution losses can be expressed, by using (2), as 

-^losses := X] \^e\'^^{Ze) = Lu. (6) 

Given a lower bound C/inin for the voltage magnitudes, we can 
therefore formulate the following optimization problem, 

vlF Lu 



mm 

qh,heC\{0} 



subject to \uh\ > C/min, V/i e C\{0} 



(7a) 

(7b) 



where voltages u are a function of the decision variables 
qh, h e C\{0}, via the implicit relation defined by the system 
of nonlinear equations (3), (4), and (5). While we are not 
considering upper bounds on the bus voltage magnitudes, in 
the form \uh\ < ?7max. it will be clear in the following that 
they could be easily incorporated with minor modifications of 
the algorithm, at the cost of a slightly more complex notation. 

In order to adopt a compact notation for the system inputs, 
measured outputs, and state, we introduce the following block 
decomposition of the vector of voltages u 



u = 



Uo 
UG 
Ul 



C"- 



^ are the 
^ are the 



where uq is the voltage at the PCC, uq S 

voltages at the microgenerators, and e 

voltages at the loads. 

Similarly, we also define = Pg+JQg and Sl = Pl+Jql- 
From a system-wide prospective, the control problem that 

we are considering is therefore characterized by 
• the input variables qg. 



the measured output variables 



Uq 
Uq 



• the unmeasured disturbances p^, qi^, pg- 

Remark. While the decision variables of the ORPF problem 
(i.e. the input variables qc) do not include the reactive power 
provided to the distribution grid by the PCC (i.e. qo = uoioA 
this quantity will also change every time the reactive power 
injections of the generators are updated by the algorithm, 
because the inherent physical behavior of the slack bus (the 



PCC) ensures that equations (3), (4) and (5) are automatically 
satisfied at every time. 

V. Dual decomposition 

In this section, in order to derive a control sttategy to solve 
the ORPF problem, we apply the tool of dual decomposition to 
(7). While problem (7) might not be convex in general, we rely 
on the results presented in [13] which show that zero duality 
gap holds for the ORPF problems, under some conditions that 
are commonly verified in practice. Based on this result, we use 
an approximate explicit solution of the nonlinear equations (3), 
(4), and (5), to derive the a dual ascent algorithm [14] that 
can be implemented by the agents. 

In order to present the approximate solution, we need the 
following technical lemma. 

Lemma 1 (Lemma 1 in [15]). Let L be the weighted Laplacian 
of Q. There exists a unique symmetric, positive semidefinite 
matrix X e R"^" such that 



XL = 7- 11^ 
Xlo = 0. 



(8) 



The matrix X depends only on the topology of the grid 
power Unes and on their impedance (compare it with the 
definition of Green matrix in [16]). 

The matrix X has some notable properties, including the 
fact tiiat 

Xhh>Xhk>0 h,keV, (9) 



and the fact that 



{ih-iky'x{ih-ik) = \zfk\, h,kGV, 



(10) 



where Z'f^ represents the effective impedance of the power 
lines between node h and k. Notice that, if the grid is radial 
(i.e. 5 is a tree) then is simply the impedance of the only 
path from node h to node k. 

By adopting the same block decomposition as before, we 
have 

"0 

M 






N 

Q 



with M e m"' 

Tij (n— m) X (n—m) 



-l)x(m-l)^ g 



-1)x(t^ 



(11) 

and e 



The following proposition provides the approximate relation 
between the grid voltages and the power injections at the 
nodes. 

Proposition 1. Consider the physical model described by the 
set of nonlinear equations (1), (2), (4), and (5). Node voltages 

then satisfy 
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Un 
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Q. 
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SG 

Sl 
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(12) 
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where the little-o notation means that lucuj^. 



o{f{UN)) _ n. 
f{UN) ~ 



Proof: The proposition descends directly from Proposi- 
tion 1 in [15]. ■ 

The quality of this approximation relies on having large 
nominal voltage Un and relatively small currents injected by 
the inverters (or supplied to the loads). This assumption is 
verified in practice, and corresponds to correct design and 
operation of power distribution networks, where indeed the 
nominal voltage is chosen sufficiently large (subject to other 
functional constraints) in order to deliver electric power to the 
loads with relatively small power losses on the power lines. In 
1 15 1, a brief discussion about how this approximation extends 
the DC power flow model [17, Chapter 3] to the lossy case, 
has been provided. 

Given the approximate exphcit expression for voltages u 
presented in Proposition 1, we proceed by specializing the 
dual-ascent algorithm to problem (7). 

The Lagrangian of the problem (7) is 



C{qG, A) = u'^Lu + (61 - vg) 



(13) 



where A is the Lagrangian multipher (i.e., the dual variables 
of the problem), where b = U^^/U^, and where vq are the 
normalized squared voltage magnitudes at the generators 



VG 



UG UG 



being the operator for element-wise multiplication. Notice 
that, in order to simplify the derivation of the algorithm, 
we squared and normalized the left and right hand side of 
the inequalities (7b), without this having any effect on the 
optimization problem. Again, ug and Vg are implicit functions 
of the decision variables qg, even if the dependence has been 
omitted in the notation. 

A dual ascent algorithm consists in the iterative execution 
of the following alternated steps 

1) minimization of the Lagrangian with respect to the 
primal variables qg 

qG{t+l) = aTgmm C{qG{t),X{t)), (14) 

2) dual gradient ascent step on the dual variables 

d£{qG{t + l),X{t)) 



\{t+l) 



dX 



(15) 



where the [•]+ operator corresponds to the projection on 
the positive orthant, and where 7 is a suitable positive 
constant. 

In order to derive an expression for the minimizer qG{t + l) 
in (14), we need the following lemma. 

Lemma 2. Consider the Lagrangian C{qG, X) defined in (13). 
The partial derivative with respect to the primal variables qG 

is 



dC{qG,X) 
dqG 



-2- {MqG + NqL - sin OMX) + o { 



Proof: From (13) we have that 

djC{qG,X) du^Lu f dvG\ 
dqG dqG \dqG J 



X. 



(16) 



In order to derive ^ g^^" , we introduce the orthogonal decom- 
position u = {u' + ju")eJ^'''+'^\ with u',u" € K". We then 
have that, via Proposition 1, 

u' = ^ («e-^(*+^)) 
and similarly 
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0" 
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PG 
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1 



= — sin^f/jvl 
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1 



By using the fact that u^Lu = u'^Lu' + u"^Lu" , we have 

dv?" Lu 



dqG 



.oqc, 





M TV 

A^^ Q 
1 



^[0 M N]L 



'N 

^{MqG+NqL)+o. 



Lu" 

qG 



1 

7/2" 



^ N 



(17) 



where we used the fact that LI = and that, by Lemma 1, 
LX = 7 - lol"^. 

The same approximate solution (12), via some algebraic 
manipulations, allows us to express vg as 

VG = \ + {e^^MsG + e^^NsL) + o . (18) 



TP 



(19) 



We therefore have that 
dvG 2 

^ = r^smOM^-o 
dqG Ufj 

and finally, from (16), (17), and (19), 

^^ = 4(M.. + iV,.-sin.MA)+o^^, 



Given Lermna 2, we propose the following update step for 
qG to minimize the Lagrangian with respect to the primal 
variables 

qa{t + 1) = qG{t) + smeX{t) - {qG{t) + M-^Nq^). (20) 

Indeed, by plugging (20) into the expression for the partial 
derivative of the Lagrangian with the respect to qG, provided 



6 



in Lemma 2, we obtain 

a£(gG(t + l),A) ^ 
dqo 

= ^ {Mqcit + 1) + Nql - smeMX{t)) + o 



TP 



= o 



r/2 



(21) 



which vanishes for large nominal voltage Un- 

In order to derive the dual ascent step for the Lagrangian 
(13), it is enough to evaluate |^ to see that 



dX 



bl - vg, 



in accordance to the well known result in dual decomposition 
that the dual ascent direction is given by the constraint 
violation. The proposed dual ascent step will therefore simply 
be 

X{t + 1) = [X{t) + 7 (W - VG{t + 1))]+ , (22) 

for some suitable 7 > 0, where we recall that vq is defined 
element- wise as Vh = \uh\'^/U'^. 

Remark. It is well known that the dual ascent method guaran- 
tees primal feasibility of the solution, but does not guarantee 
primal feasibility of each intermediate step. Therefore, in the 
scenario that we are considering, voltage constraints (7b) will 
not be satisfied at every time. This is however acceptable: volt- 
age constraints can be intended as soft constraints, because 
they do not derive from physical constraints on the system (like 
Kirchhoff laws). It will be shown via simulations in Section IX 
how the constraint violation is actually quite limited, also in 
the case of time-varying loads. 

VI. Synchronous algorithm 

In this section, we show how the agents can implement the 
primal and dual steps (20) and (22) derived in the previous 
section. We assume here that the agents are coordinated, i.e. 
they can update their state variables and Xh, h e C\{0}, 
synchronously. 

In order to derive the update law for the agents, we need to 
introduce the following matrix G. 

Lemma 3. There exists a unique symmetric matrix G G 
R^'X"' such that 




M 
Gl = 0. 



G = I 



11? 



Proof: The following synmietric matrix G satisfies the 
conditions. 



G 



-M-^l 



(23) 



The proof of uniqueness, that we omit here, follows exactly 
the same steps as in the proof of Lemma 1. ■ 



The matrix G has also a remarkable sparsity pattern, as the 
following lemma states. 

Lemma 4. The matrix G has the sparsity pattern induced by 
the Definition 1 of neighbor agents in the cyber layer, i.e. 



keAf{h). 



The proof is provided in the Appendix A, where we also 
discuss how the elements of G can be estimated by the agents, 
given a local knowledge of the power grid topology and 
parameters. 

We therefore propose the following synchronous algorithm. 

Synchronous algorithm 

Let all agents (except the PCC) store an auxiliary scalar 
variable A^. Let 7 be a positive scalar parameter, and let 9 
be the impedance angle defined in Assumption 1. Let Ghk be 
the elements of the matrix G defined in Lemma 3. At every 
synchronous iteration of the algorithm, each agent h € C\{0} 
executes the following operations in order: 

• gathers the voltage measurements 

{uk = \uk\exp{jZuk), k e J\f{h)} 

from its neighbors; 

• updates the auxiUary variable Xh as 



Xh + 1 



^ N 



^ N 



(24) 



• based on the new value of Xh, updates the injected 
reactive power qh as 

qh ^ qh + sin eXh 

+ X] Ghk\uh\\uk\ sm{Zuk - Z.Uh - 0). 

(25) 

It can be shown, by using Lemma 4 and via some algebraic 
manipulation, that the update (25) can be also rewritten as 

qc ^ qG{t)+sm0X{t) + Q (^e"^'^ [O diag(uG)] ^° ) > 

which, by using the expression for u provided by Proposition 
1, is equal to 



qG<-qG + sineX -{qc + M ^NqL) + o 



1 



Similarly to what has been done in (21), we have that after 
the update 

djC{qG,X) 



dqo 



1 

IP 



and therefore the update minimized the Lagrangian with 
respect to the primal variables, up to a term that vanishes 

for large Un- 

It is important to notice that the proposed synchronous 
algorithm requires that the agents actuate the system at every 
iteration, by updating the amount of reactive power that the 
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agent command to the microgenerators. Only by doing so, the 
subsequent measurement of the voltages will be informative 
of the new state of the system. The resulting control strategy 
is thus a feedback strategy that necessarily requires the real- 
time interaction of the controller (the cyber layer) with the 
plant (the physical layer), as depicted in Figure 4. This tight 
interaction between the cyber layer and the physical layer is 
the fundamental feature of the proposed approach, and allows 
to drive the system towards the optimal configuration, which 
depends also on the reactive power demand of the loads, 
without collecting this information from them. In a sense, 
the algorithm is inferring this hidden information from the 
measurement performed on the system during its execution. 

VII. Asynchronous algorithm 

In order to avoid the burden of coordination among the 
agents, we also propose an asynchronous version of the 
algorithm, in which the agents corresponding to the micro- 
generators update their state {q^, A^) independently one from 
the other, based on the information that they can measure and 
that they can gather from their neighbors. 

We assume that each agent (except for the agent located at 
the PCC) is provided with an individual timer, by which it is 
triggered. Timers tick randomly, with exponentially, identically 
distributed waiting times. 

Asynchronous algorithm 

Let all agents (except the PCC) store an auxiliary scalar 
variable Xh- Let 7 be a positive scalar parameter, and let 9 be 
the impedance angle defined in Assumption 1 . Let Guk be the 
elements of the matrix G defined in Lemma 3. 

When agent h e C\{0} is triggered by its own timer, it 
performs the following actions in order: 

• gathers the voltage measurements {uk = 
Iwfcl exp(jZ?ife), k E Af{h)} from its neighbors; 

• updates the auxiliary variable as 

12 ~ 



Ah + 7 



C/2. 

nun 
2 
JV 



U 



\uh\ 
TP 



• based on the new value of A/,, updates the injected 
reactive power as 

Qh qh + sin OXh 

+ ^ Ghk\uh\\uk\sm{Zuk - ^Uh - 0). 

The update equations for the asynchronous algorithm are 
exactiy the same of the synchronous case. Here, however, we 
update both the primal and the dual variable of the agents 
independently and asynchronously. 

VIII. Convergence analysis 

In this section, we study the convergence of both the 
synchronous algorithm proposed in Section VI and of the 
asynchronous algorithm proposed in Section VII. For the 
analysis of the stability of both of them, we adopt again the 
approximated model proposed in Proposition 1, and we neglect 
the infinitesimal terms. 



Notice that, by neglecting the infinitesimal terms, both the 
voltages u and the squared voltage magnitudes vq become 
affine functions of the decision variables qc, and therefore the 
ORPF problem (7) becomes a convex quadratic problem with 
linear inequality constraints, for which strong duality holds. 

Let us define 

WG := 1 + [ei'^MsG + c^'^Nsl) , (26) 



U 



N 

For the synchronous case we consider the update equations 

Xit + l) = [X{t)+j{bl-VGm+, (27) 

for the dual variables, and 

qait + l) = qGit)+sm0X{t+l)-{qG{t) + M-'NqL). (28) 

for the primal variables. 

Notice that the equilibrium {q^, A*) of (27)-(28) is charac- 
terized by 

bl-VGit)<0 and sindX* - {qQ + M-'^NqL) = 0, 
and therefore, via (18) and via Lemma 2, we have that 



and 



dC{q*G,X*) _ fj_ 
dq*G ''[Ul I' 



' N 



which correspond, up to infinitesimal terms, to the necessary 
conditions for the optimality of (g^, A*) according to Uzawa's 
saddle point theorem [18]. 
The following convergence result holds. 

Theorem 1 (Synchronous case). Consider the dynamic system 
described by the update equations (27) and (28). The equilib- 
rium {q*,X*) is asymptotically stable if 



7< . o 



TP 



sm^6{m-l)D' 

where m — 1 is the number of microgenerators and D := 
max/j \Z^\ is the maximum electric distance of a microgen- 
erator from the PCC. 

The proof is presented in Appendix B. 
For the asynchronous case, on the other hand, we introduce 
the following assumption. 

Assumption 2. Let {T^^^}, i e N, be the time instants in 
which the agent h is triggered by its own timer. We assume that 
the timer ticks with exponentially distributed waiting times, 
identically distributed for all the agents in C\{0}. 

Let us define the random sequence h{t) G C\{0} which tells 
which agent has been triggered at iteration t of the algorithm. 
Because of Assumption 2, the random process h{t) is an i.i.d. 
uniform process on the alphabet C\{0}. We therefore consider 
the following update equations, in which only the component 
h{t) of the vectors qG and A is updated at time t: 



hinit + 1) = [^hit){t) + 7il(t) (w - Mt)) 

Xk{t +1) = Xk{t + 1) for all A; ^ h{t), 



(29) 
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Figure 4. A block diagram representation of tlie synchronous control algorithm proposed in Section VI, where the tight interconnection of the cyber and the 
physical layer (i.e. the feedback strategy) is evident. 



and 

qh{t) (i + 1) = Qh{t) {t) + sin 9Xh(t) + 1) 

-l^(,)((?G(i) + M-iiVgi), (30) 
Qkit + 1) = qk{t) fovidl k^h{t). 

Notice that, also in the asynchronous case, Uzawa's neces- 
sary conditions for optimality are satisfied at the equilibrium 

of (29)-(30). 

The following convergence result holds. 

Theorem 2 (Asynchronous case). Consider the dynamic sys- 
tem described by the update equations (29) and (30). Let 
Assumption 2 hold. Then the evolution t — > {q{t),\(t)) 
converges almost surely to the equilibrium (g*, A*) if 



~ sm^6{m-l)D' 

where m — 1 is the number of microgenerators and D := 
max/i l-^oll maximum electric distance of a microgen- 

erator from the PCC. 

The proof is presented in Appendix B. 

IX. Simulations 

The algorithm has been tested on the testbed IEEE 37 
[8], which is an actual portion of 4.8kV power distribution 
network located in California. The load buses are a blend 
of constant-power, constant-current, and constant-impedance 
loads, with a total power demand of almost 2 MW of active 
power and 1 MVAR of reactive power (see [8] for the testbed 
data). The length of the power lines range from a minimum 
of 25 meters to a maximum of almost 600 meters. The 
impedance of the power lines differs from edge to edge (for 
example, resistance ranges from 0.182 Vllkm to 1.305 0/km). 
However, the inductance/resistance ratio exhibits a smaller 
variation, ranging from Az^ = 0.47 to /-Z^ = 0.59. This 
justifies Assumption 1, in which we claimed that Zzg can 



PCC 




Figure 5. Schematic representation of the IEEE 37 test feeder [8], where 5 
microgenerators have been deployed. 



be considered constant across the network. We considered the 
scenario in which 5 microgenerators have been deployed in 
this portion of the power distribution grid (see Figure 5). 

The lower bound for voltage magnitudes has been set to 
96% of the nominal voltage. Both the synchronous and the 
asynchronous algorithm presented in Section VI and VII have 
been simulated on a nonlinear exact solver of the grid [19]. 
The approximate model presented in Proposition 1 has not 
been used in these simulations, being only a tool for the 
design of the algorithm and for the study of the algorithm's 
convergence. The parameter 7 has been chosen as one half 
of the bound indicated by Theorem 1 and Theorem 2 for 
convergence. Timers in the asynchronous case have been tuned 
so that each agent is triggered, in expectation, at the same rate 
of the synchronous case. 
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A time-varying profile for the loads has been generated, in 
order to simulate the effect of slowly varying loads (e.g. the 
aggregate demand of a residential neighborhood), fast chang- 
ing demands (e.g. some industrial loads), and intermittent large 
loads (e.g. heating). 

The results of the simulation have been plotted in Figure 
6. For both the synchronous and the asynchronous case, the 
power distribution losses and the lowest voltage magnitude 
measured by the microgenerators, are reported. The dashed 
line represents the case in which no reactive power compen- 
sation is performed. The thick black line represents the best 
possible strategy that solves the ORPF problem (7) (computed 
via a numerical centralized solver that have access to aU the 
grid parameters and load data). The thin red line represents 
the behavior of the proposed algorithm. 

It can be seen that the proposed algorithm achieves prac- 
tically the same performance of the centralized solver, in 
terms of power distribution losses. Notice however that the 
proposed algorithm does not have access to the demands of the 
loads, which are unmonitored. The agents, located only at the 
microgenerators, can only access their voltage measurements 
and share them with their neighbors. Notice moreover that, 
as expected for duality based methods, the voltage constraint 
is satisfied only at the steady state. Therefore, in the time 
varying case simulated in this example, the voltage sometimes 
falls lightly below the prescribed threshold, when the power 
demand of the loads present abrupt changes. This effect is 
even more noticeable in the asynchronous case, where iteration 
times are not evenly spaced. It should be remarked, however, 
that the extent of this constraint violation depends on the 
rate at which the algorithms are executed, compared with the 
rate of variation of loads, and is ultimately a function of the 
communication resources that are available at the cyber level. 

X. Conclusions 

In this paper we proposed a distributed control law for 
optimal reactive power flow in a smart power distribution grid, 
based on a feedback strategy. Such a strategy requires the 
interleaving of actuation and sensing, and therefore the control 
action (the reactive power injections qh,h £ C\{0}) is a func- 
tion of the real time measurements (the voltages Uh,h G C). 
According to this interpretation, the active power injections 
in the grid {ph-,h G V) and the reactive power injection of 
the loads {qh,h G V\C) can be considered as disturbances 
for the control system. As it happens in all feedback control 
systems, these quantities do not need to be known to the 
controller: in some sense, the agents are implicitly inferring 
this information from the voltage measurements performed on 
the grid. Similarly, it is also well known that the presence 
of feedback in the control action makes the closed loop 
behavior of the system less sensitive to model uncertainties. 
While we have not elaborated this point, we expect that 
the proposed strategy exhibits this kind of robustness. These 
features differentiate the proposed algorithm from most of 
the ORPF algorithms available in the power system literature, 
with the exception of some works, like [20], where however 
the feedback is only local, with no communication between 



the agents, and of [6] and [7]. Moreover, in the proposed 
feedback strategy, the controller (or optimizer) does not need 
to solve any model of the grid in order to find the optimal 
solution. While some knowledge of the grid parameters are 
used in the design and tuning of the algorithm, the online 
controller does not need to know or solve the nonlinear grid 
equations that are generally a critical issue in offline ORPF 
solvers. On the contrary, the computational effort required for 
the execution of the proposed algorithm is minimal. These 
features are extremely interesting for the scenario of low 
voltage or medium voltage power distribution networks, where 
real time measurement of the loads is usually not available, 
the grid parameters are partiaUy imknown, and many buses are 
unmonitored. 

While a feedback approach to the ORPF problem is quite 
a new contribution, similar methodologies have been used to 
solve other tasks in the operation of power grids (see Figure 
7). In particular, in order to achieve realtime matching of the 
power demand and power generators, synchronous generators 
are generaUy provided with a local feedback control that 
adjusts the input mechanical power according to frequency 
measurements, in order to keep frequency sufficiently close 
to the desired nominal value (the primary frequency control 
in [21], see also [22], [23]). By adding a communication 
channel (a cyber layer) that enables coordination among the 
different agents, it is also possible to drive the system to 
the configuration of minimum generation costs [24]. Notice 
that, in this scenario, generators do not have access to the 
aggregate active power demand of the loads, but infer it 
from the purely local frequency measurements. This example 
share some qualitative similarities with the original approach 
presented in this paper, where agents do not have access to 
the individual or aggregated reactive power demand of the 
loads, but however can use voltage measurements to infer this 
information and, by coordinating their action, to drive the grid 
to the configuration of minimum losses. 

As suggested in [25], a control-theoretic approach to the 
problem of optimal power flow enables also a number of 
analyses on the performance of the closed loop system that 
are generally overlooked when tackling the problem with 
the tools of nonUnear optimization. Examples are L2-like 
metrics for the resulting losses in a time-varying scenario 
(see for example the preliminary results in [26 1), robustness 
to measurement noise and parametric uncertainty, stability 
margin against communication delays. These analyses, still 
not investigated, are also of interest for the design of the 
cyber architecture that has to support this and other real time 
algorithms, because they can provide specifications for the 
communication channels, the communication protocols, and 
the computational resources that need to be deployed in a 
smart distribution grid. 

Appendix A 

G-PARAMETERS 
Let us define the following parameters ghk, h,k G C. 
Deflnition 2 (^-parameters). For each pair h,k G C, let us 
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Figure 6. For both the synchronous and the asynchronous case, power distribution losses and the lowest measured voltages are plotted for the following 
cases: when no reactive power compensation is performed (dashed line), when an ideal centralized numerical controller commands the microgenerators (thick 
black line), and for the proposed algorithm, where microgenerators are commanded via a feedback law from the voltage measurements (thin red line). 
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Figure 7. Qualitative comparison between primary frequency control in 
power grids (first panel) and the proposed feedback strategy for optimal 
reactive power flow (second panel). 



define the parameter 

9hk = \ih\ 



Uk 
Ui 

k 



= 1 

= 0, £ e C\{k} 
^Q, lie 



(31) 



i.e. the current that would be injected at node h if 

• node k was replaced with a unitary voltage generator; 

• all other nodes in C (the other agents) were replaced by 
short circuits; 

• all nodes not in C (load buses) were replaced by open 
circuits. 

Notice that the parameters g/j^. depend only on the grid 
electric topology, and that 

9hk 7^ if and only if fee JV{h). 

Figure 8 gives a representation of this definition. Notice 
that, in the special case in which the paths from h to its 



ui—O 

yeec\{h} 



k e U(h) 



Uh = 1 



k e N(h) 
<9 — 



U£ — 

V£ e c\{h} 



lift =1 
— ©— 



il=0 











Figure 8. A representation of how the elements G^h are defined. Notice 
that in the configuration of the left panel, as the paths from h to its neighbors 
k £ Af{h) do not share any edge, the gains Gf^^ corresponds to the absolute 
value of the path admittances 



neighbors are all disjoint and unique paths, then Ghk ~ 
(SeGPhfc Ne|)~^, i-e. the inverse of the impedance of the path 
connecting h to k. 

As suggested in [12], these parameters can be estimated in 
an initialization phase via some ranging technologies over the 
PLC channel. Alternatively, this limited amount of knowledge 
of the grid topology can be stored in the agents at the 
deployment time. Finally, the same kind of information can 
be also inferred by specializing the procedures that use the 
extended capabilities of the generator power inverters for 
online grid sensing and impedance estimation [27], [28]. 

The following Lemma shows that the elements g^k in 
Definition 2 correspond to the elements Ghk of the matrix 
G in Lemma 3. 

Lemma 5. Let G be the matrix defined in Lemma 3. Then, 
for all h, k, Ghk = ghk o.s defined in Definition 2. 

Proof Let G be the matrix whose elements are the 
parameters ghk- From Definition 2, we have that when — 0, 



= cxp(-j0)G 



"0 
UG 



(32) 



From circuit theory considerations, this implies that Gl = 0. 
From (3) and by using the matrix X defined in Lemma 1, if 
= 0, we have 
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UG 





M 



(33) 



Therefore, by plugging (32) into (33), we have that 
{l-lll) = [li,]G. m 



Appendix B 
Proof of Theorems 1 and 2 

We first introduce the following technical lemma. 

Lemma 6. Let M be the (m — 1) x (m — 1) block of X as 
defined in (11). Let p{M) be its spectral radius. Then 

p{M) <{m- l)D 

where D max/j \ Z'^Qh\ ''^ maximum electric distance of 
a microgenerator from the PCC. 
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Proof: Because of (9), we have that M^h > Mhk > 
0, V/i, k. Therefore 



p{M) 



max 



hkVk 



< max 
ll'/INi 



Ef^'^'^E^l 

h \ k / 



^ JE (MhhVn^) 



< 



^ E ^maxM/j^Vm - 1 



maxMft,j(m - 1). 
ft 



|Zgf |, we prove the 



By using the fact that, via (10), M^h 
statement. ■ 
Proof of Theorem 1: Let {qq, A*) be the equilibrium of 
the system (27)-(28), which satisfies the following equations 



sin^A* 



a;; = 0.. 



(34a) 
(34b) 
(34c) 



{q*G + M-'NqL)^0 

b - Vh{q*) < V/i e C 
b - Vh{q*) < <^ 

We introduce the two following quantities 

x{t) = qcit) - qh and y{t) = \{t) - X* . 
The update for x{t) is, given (28), 

x{t + 1) ^ qait + 1) - q* 

= x{t) + sin^A(t + 1) - {qcit) + M-^NqL) 

-smOX* +q*G + M-^NqL 
= sm.ey{t+l), 

where in the second step we used condition (34a). 

We then consider the update equation for y{t). From (27) 
we have 

y{t + i) = 

= X{t + 1) - A* 

= [A(t) + 7 (61 -%(*)))]+ -[A1 

y{t) + A* + 7 (61 - v{q*)) - 7-^ smeMx{t) 



-[A1+, 
where we used the fact that 

v{q{t)) = v{q*) 



U 



N 



2 smeMx{t). 

N 



Let us define the constant vector 
a = A* +7 (61 
We then have, by using (34b) and (34c), 



2/(i+l) 



y{t) + a - 7—2- sm0Mx{t) 



which, by using the update equation for x, becomes 

2 



y{t + 1) 



/ - 7-^ sin^ 0M ) y{t) + a 



- \a\. 



By using the fact that 

||a+-6+||<||a-6||, 

we have that 

2 



\\y{t + m< 



and therefore, as M is symmetric and positive definite, y{t) 
converges to zero if 



r/2 



^ ^ sin^ ep{M) ' 

where p(M) is the spectral radius of M. 

By using Lemma 6, we then finally have the sufficient 
condition for convergence 



7< - 



r/2 



sin^ 6»(rn- l)maxft ' 



Proof of Theorem 2: Consider the update equations (29) 
and (30) for the dual variables A and for the primal variables 
qc, respectively. Introducing the same variables y and x as in 
the proof of Theorem 1, we can write, assuming, without loss 
of generality, that node h is the node performing the update 
at the t-th iteration 

yh{t + i) = 

= Xh{t + l)-Xl 

= [Xn(t)+^ll{bl-VG{t))\^-[Xl]^ 

= [yh{t) + XI+ 7lft (61 - iait))]^ - [Xl] + 

= [yh{t) + xi + ^{bh- Mq*)) - 7 {vhit) - vh{<i*))]+ 
- [Kh 

Let ah = XI + "f {bfi - Vh{q*)) and observe that [A^]_(_ = 
[aft,]_,_. Hence we can write 



yh{t + i) = 

= [yh{t) + a/i - 7 {vh{t) - Vh{q*))]+ - M+ 

= yh{t) + ah-l^smeilM{q{t)-q 



yh{t) +ah- 7772- sin Oli Mx{t) 



[ah 



where in the second equality we have used the fact that Vh{t) — 
Vh{<t) = ^s^rieilM{q{t)-q*). 
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As far as the variable x is concerned we have that 

Xh{t+l) = 

= qh{t+l)-ql 

= qh(t) + sin0Aft(t + 1) - ll{qG{t) + M'^NqL) - ql 
= qh{t) -ql + sindXhit + 1) 

- lliqcit) -qh + qh + M-^NqL) 
= Xh{t) + smeXh{t + I) - llx{t) 

-ll{q*G + M-'NqL) 
= sineXhit + 1) - llM-\Mq*a + TV^l) 

Now from (34a) we have that sin 6X1 = '^h {ih + M-'^NqL) 
and, in turn, 

Xh{t+l) = sinOXhit + 1) - sin^A^ 
= smeyh{t + l) 

Observe that Assumption 2 implies there exists almost surely 
a positive integer T such that any node has performed an 
update within the window [0,T]. It follows that, for t > T, 
x{t) = sin 9y{t). Accordingly for t > T the update for y can 
be rewritten as 



yh{t + i) 



yh{t) +Oih- 7772- sin" OllMy{t) 



• TJ2 



and, clearly, yk{t + 1) = jyfc(i), for k ^ h. 

Now let I - 7^ sin^ Olhl^M. Using the fact that 

||a+ - < \\a - b\\ it follows 

ll?y(t + i)ll < WPkvm- 

Consider the evolution of the quantity E From 
the above inequahty we get that 

E[||y(t + l)f] <E[y{tfP^P,y{t)] 

= traceE [y{tf P^ Pnyit)] 

= trace {E [PjP^y(%(if ] } 

< ||E[PjP;.]|| trace {E[y(%(tf]} 

<||E[P,^P,]||E[||y(i)ir] (35) 

where, given two semidefinite matrices A, B, we used the fact 
that trace {AB} < || A||trace {B}. Let us compute E [P/f P?i] . 
Observe that 



PhPh = 



IP 



= 1-1^ sin^ 0{Mlhli + Ih^M) 

'J N 

from which we get, by using Assumption 2, 



m 
I 



4 sin^ e 2 4 sin'' , 

m — 1 U% m — 1 



By straightforward calculations one can see that the spectral 
radius of the above matrix is smaller than 1 if 



7< 



r/2 



sin^ ep{M) 



where we recall that p{M) denotes the spectral radius of M. 
Then, by Lemma 6, we have 



7 < 



TP 



sin^eim- l)maxh\Z^l\' 



(36) 



The proof can be concluded by invoking the Supermartin- 
gale Convergence Theorem [29]. For the sake of the clarity, 
we recall that this theorem states that, if X^, t > 0, is a 
nonnegative random variable such that E[Xi] < +oo and if 
E[Xt+i|J^t] < Xt with probability one, where J^t denotes the 
history of the process Xt up to time t, then Xt tends to a limit 
X with probability one, and limt^ooE[Xt] = E[X]. 

Let Xt = Observe that 

E [\\y{t + l)f I y{t)] < E [y{tf P/^ P,Mt) 1 2/ W] 



< 



< 



l|E [PhPh] I 



urn' 



where the last inequality follows from the fact that 
||E [P/f P/i] II < 1 under condition (36). Moreover in- 
equality (35) together with ||E [PjP/i] || < 1 implies that 
limt^ooE = 0. Hence E converges with 

probability one to a nonnegative random variable X such that 
£ [X] = 0, i.e., X is the nuU random variable. This concludes 
the proof. ■ 
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